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1 Introduction 


During this contract period, our work has focused on improvements to elliptic 
grid generation methods. There are two principle objectives in this project. 
One objective is to make the elliptic methods more reliable and efficient, and 
the other is to construct a modular code that can be incorporated into the 
National Grid Project (NGP), or any other grid generation code. Progress 
has been made in meeting both of these objectives. The two objectives are 
actually complementary. As the code development for the NGP progresses, 
we see many areas where improvements in algorithms can be made. 


2 Current Research 

During the past year, most of our research has been on the development of the 
elliptic method for generating surface grids. This has included the derivation 
of the orthogonality and control function options that are normally used with 
two and three-dimensional elliptic methods. These options were not available 
in the EAGLE grid generation code and, as far as is known, had not been 
discussed in any previous papers on elliptic surface grid generation. The 
convergence problems which were initially encountered with the surface code 
have been reduced by using a new form of upwind differencing on the elliptic 
equations and by reducing the size of the NURBS control net when computing 
the surface derivatives. However, convergence can still be a problem if the 
surface is poorly parameterized. Much of the success in the development of 
surface grid generation methodology has been due to the use of the standard 
NURBS representation of surfaces. The present state of development will be 
present at the International Conference on Hydro-Science and Engineering in 
Washington, DC. A copy of the paper ” Elliptic Grid Generation on Surfaces”, 
which will appear in the Proceedings, has been attached to this report. 

The major task in the development of the elliptic volume code has been 
removing the dependence on the EAGLE data structure and writing a driver 
routine that can accept and process arbitrary grid blocks. At the same time 
we have improved some of the algorithms dealing with control functions and 
boundary orthogonality. There were several places where iterative procedures 
could be replaced by direct computations without a noticeable decrease in 
quality of the resulting grid. 



The main elliptic volume and surface grid generation subroutines have 
been incorporated into the NGP system and are presently operational. The 
subroutines for the various control function and orthogonality options are also 
in the system, but they cannot be used since the panels for picking the options 
have not been built into the graphical user interface. However, the basic 
structure of the subroutines have been finalized so that code development 
can continue. 



To be presented at the International Conference on Hydro-Science and 
Engineering, Washington, DC, June 1993. 
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ABSTRACT 

Elliptic grid generation methods for the construction of surface grids are discussed. New 
techniques for computing the control functions that determine the grid point distribution 
and for imposing boundary orthogonality are developed. Sample grids are presented to 
demonstrate the application of these methods in the solution of some common problems in 
grid generation. 

INTRODUCTION 

Grid generation is a critical area for many numerical simulation problems. Whether 
using structured or unstructured grids, it is up to the numerical analyst to distribute points 
throughout an often complex physical region. The distribution of the grid points and t e 
geometric properties of the grid such as skewness, smoothness, and cell aspect ratios have a 

major impact on the accuracy of the simulation. 

Most grid generation proceeds in stages. In two dimensions, the grid is first constructed on 
the boundary curve or curves and then constructed on the interior of the physical region. In 
three dimensions, the grid generation process proceeds from curves to surfaces and then to the 
interior of the physical region. At each stage of the grid generation process, the construction 
of the desired grid often follows in two steps. First, a grid is constructed by interpolation 
from the boundary of the region or surface, and then this grid is smoothed, and possibly 
modified in other ways, by an iterative procedure. The most successful iterative smoothing 
schemes are based on elliptic systems of partial differential equations that relate the physica 
and computational variables. The elliptic system may be applied to the boundary grids, the 
interior grids, or both. The elliptic system may preserve the original distribution of grid 
points or redistribute points as when constructing an adaptive grid. Orthogonality of the 
grid may be imposed along certain boundary components of the physical region. 

The most difficult and less developed area of elliptic grid generation is that of grid 
generation on surfaces. This is because the grid must not only be smoothed, but the grid 
points must also stay on the surface. For surfaces defined by parametric equations, the 
simplest technique for achieving this goal is to work in parameter space rather than in 
the physical variables of the surface. There are some disadvantages associated with this 
approach. The differential equations become more complicated and contain two sets of 
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derivatives, the derivatives of the physical variables with respect to the parametric variables 
and the derivatives of the parametric variables with respect to the computational variables. 

This report will discuss some of the current problems in elliptic grid generation on surfaces 
and describe some new features and enhancements that have been recently developed. Grid 
orthogonality has been one area that has been of particular interest. Although orthogonal 
grids are generally not necessary for most numerical algorithms, there must be some limit on 
the skewness of the grid. Boundary orthogonality is also desireable whenever it is necessary 
to impose Neumann type boundary conditions on one or more of the physical variables in the 
problem. Also, some of the popular ways of treating Neumann boundary conditions are not 
accurate when the grid is not orthogonal and others become unstable if the grid is extremely 
skewed near the boundary. 

ELLIPTIC SURFACE EQUATIONS 

The elliptic equations for surface grids are a generalization of Laplace equations for 
harmonic mappings of plane regions. Let S be a simply-connected surface defined by the 
parametric equations 

T — f(u,v), 0<ti,v<l. 

Further, let u and v be functions of the computational variables £ and tj where 0 <£,?/ < 1. 
Now a cartesian coordinate system in the computational rectangle generates a curvilinear 
coordinate system in the parametric rectangle which maps to a curvilinear coordinate system 
on the surface. Thus a uniform grid in the computational rectangle generates a curvilinear 
grid on the surface. The elliptic system of equations which relates the parametric and 
computational variables is given as 
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Although we will not attempt to derive these equations, they are related to conformal map- 
pings onto surfaces and a discussion may be found in Mastin and Thompson (1984). Also, 
some of the theoretical results from harmonic mappings of plane regions can be verified for 
solutions of the above system. These same equations may also be derived from a differential 
geometric point of view as in the paper by Warsi (1990). The functions P and Q are used to 
control the grid point distribution and would not appear in the derivation from conformal 
or harmonic mappings. 



The system of partial differential equations may be solved with either Dirichlet or Neu- 
mann type boundary conditions depending on whether the grid points on the boundary of 
the surface are to remain fixed or are allowed to move during the construction of the grid. 
Imposing boundary orthogonality is the most common case where grid points are allowed to 
move along the boundary. However it is orthogonality of the grid on the surface and not in 
the parametric region that is desired. Thus the appropriate boundary conditions must be 
imposed on the elliptic system in the parametric region so that the grid is orthogonaljm the 
surface. Two grid lines will be orthogonal on the surface if the inner product • r„ = 0. 
This can be expanded using the chain rule to give the following equation 

9n u i u n "b 9n v i v n + 9ni u i v 'i "b u v v t) ~ ^ 

The orthogonality condition can be formulated as derivative boundary conditions for the 
above elliptic system. If the boundary segments u = 0 and u = 1 are considered, then the 
orthogonality condition reduces to 

9?2 v Z "b 9\2 u i = ^ ' 

Similarly, for the segments v = 0 and v = 1, orthogonality is imposed by the equation 

9n u n + 9\2 v n = 

Thus the elliptic system together with the given values of u and v on the boundary and 
the orthogonality boundary conditions form a mixed boundary value problem. Equation (4) 
determines the values of v along the boundary where u = 0 or u = 1, and Equation (5) 

determines the values of u where i> = 0 or v = 1. 

The control functions P and Q must be selected so that the grid has the required dis- 
tribution of grid points in the computational field. Taking P = Q — 0 tends to generate 
a grid with a uniform spacing which is seldom desired. In most cases there is a need to 
concentrate points in some area of the surface such as along certain boundary contours. 
There are two methods of computing the control functions. The first is to compute P and 
Q from some initial grid. All the derivative terms in the elliptic system can be computed 
from given grid point values leaving two unknowns, P and Q , which can be determined from 
the two equations. Now these control function values are smoothed so that the final ellip- 
tic grid will be smoother and generally more orthogonal than the initial grid. The control 
functions can be computed from an initial grid only if the Jacobian of the transformation 
from computational to parametric variables is nonvanishing. Since this may not always be 
the case for interpolated grids, there is also the option of computing the control functions 
from the boundary distribution. The appropriate values for the control functions on the 
boundary can be derived by assuming the grid lines are orthogonal at the boundary and the 
spacing normal to the boundary is uniform. The development basically follows that for two 
dimensional plane regions as found in the book by Thompson et. al. (1985). The formulas 
for P and Q in the surface Equations (1) and (2) are given below. 

P = + ( 6 ) 

s i 

Q = + ( 7 ) 
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The variable s is the arc length parameter along the boundary of the surface and the variables 
1C i and JC 2 denote the curvature of the boundary curves f = constant and rj = constant , 



respectively. The curvatures are given by 



where the Christoffel symbols T £ in the previous formulas are given by 
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and 
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Not all of the terms in Equations (6) and (7) for P and Q can be evaluated on the boundary. 
It is clear that the arc length derivatives for P can be evaluated on the rj = 0 and tj = 1 
boundary curves and interpolated in the interior of the computational region, while the arc 
length derivatives for Q can be evaluated on £ = 0 and £ = 1 and interpolated in the interior. 
In a similar manner, the curvature in each control function can be computed on opposite 
boundaries and interpolated. Note that the curvatures cannot be computed solely from 
information on the boundary curves. In this manner all of the arc length derivatives and 
curvatures, and therefore P and Q, can be computed throughout the computational region. 
This second control option can be used to determine the distribution of grid points on a 
surface even when an initial grid is of such poor quality in the interior of the surface that it 
cannot be used itself to generate appropriate control functions. 

In the above discussion it was shown how derivative boundary conditions can be used to 
impose orthogonality. In a numerical algorithm, this requires moving the grid points along 
the boundary of the surface. Orthogonality can also be imposed by adjusting the control 
functions near the boundary and keeping the boundary points themselves fixed. Since there 
are two control functions, not only can the angle at the boundary be chosen, but the distance 
to the first grid line off of the boundary curve can also be chosen. If it assumed that the 
grid lines are orthogonal, the values of P and Q can be determined from the elliptic system 
(1) and (2) yielding the expressions 
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Along a boundary component where a specified grid spacing is desired, that distance is used 
for the value of g n or g 22 - Some of the other derivatives in the equations for P and Q can 
be computed from the fixed values of u and v on the boundary. Those derivatives that must 
be computed using interior values are updated during the iterative solution of the elliptic 
system. At each iteration, these boundary control functions are smoothly blended with the 
interior control functions to generate a smooth grid that is orthogonal and has a specified 
spacing along assigned boundary components. 

EXAMPLES 

Two computational examples will be presented to demonstrate the utility of the method 
in grid generation. In both cases the surface is defined as a NURBS (NonUmform Rational 
B-Spline) surface. The elliptic system is solved using a finite difference discretization and an 
SOR iterative method. In some cases, the SOR method does not converge. In those cases, 
an acceptable grid can often be generated by setting the right hand sides of Equations (1) 
and (2) to zero. Of course, curvature information about the surface is lost since all of the 
second order derivatives with respect to the parametric variables appear on the right hand 
side of the equations. In both of the examples the grids are very coarse so that the grid 

properties can be clearly seen from the plots. 

The first example is a simple surface with a complex parameterization. Two grids are 
shown in Figure 1. One grid is constructed by transfinite interpolation in the parametric 
region. The other grid is generated by smoothing the first grid with the elliptic system. 
The second grid is certainly smoother and more uniform. The control functions were taken 
to be zero and the boundary points were allowed to move so that the grid is orthogonal at 
the boundary. Some nonuniformity is noted in the grid, but these are due to small changes 
in the curvature of the surface and not to the parameterization. At this point, it is not 
known if these effects are inherent in the elliptic system or are due to truncation error in the 
discretization. 

The second example is a more realistic configuration. A space shuttle type geometry is 
depicted with an interpolated and a smoothed grid in Figure 2. Due to the large curvature of 
the surface, convergence problems were encountered in the numerical solution of the elliptic 
system. Therefore the right hand sides of the elliptic system in Equations (1) and (2) were 
set to zero. Although there is not a great deal of difference in the two grids, the elliptic 
grid does result in a more uniform distribution especially in front of the wing where the grid 
constructed by interpolation has a sparser distribution of points. 

CONCLUSIONS AND RECOMMENDATIONS 

The elliptic equations for surface grids have been presented along with the control func- 
tions and orthogonality techniques which allow a great deal of freedom in designing high 
quality grids. This capability would be most useful when the parameterization of the sur- 
face is such that interpolation of parametric values does not give a satisfactory grid because 
of a poorly parameterized surface. Since this situation arises frequently when surfaces are 
defined by CAD packages, the capability to smooth and improve surface grids is essential 
in any state-of-the-art grid generation code. The method described in this report works 
well on surfaces which have a smooth parametric representation with moderate values for 
the .parametric derivatives. If there are large variations in the parametric derivatives, then 
the elliptic system becomes difficult to solve and the standard iterative methods, such as 
SOR, will not converge. However, this is precisely the case when the elliptic grid generation 



Figure 1: Initial Grid (left) and Elliptic Grid (right) for a Surface 



Figure 2: Initial Grid (left) and Elliptic Grid (right) for a Shuttle Configuration 


methods can be most useful. Thus there is a need for more study on algorithms to solve 
nonlinear elliptic systems such as those encountered in this report. 
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